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WKB analysis of edge states in graphene in a strong magnetic field. 
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We investigate the fine structure of the edge states energy spectrum for zigzag and armchair 
ribbons of graphene in a strong magnetic field. At low energy, the spectra can be described by 
an effective Schrodinger Hamiltonian with a double well potential, symmetric in the zigzag case 
and asymmetric in the armchair case. We develop a semiclassical formalism based on the WKB 
approximation to calculate analytically the energy spectrum for the two types of edges, including 
regions which were not studied earlier. Our results are in very good quantitative agreement with 
numerical calculations. This approach leads to a qualitative description of the spectra in terms of 
the quantization of unusual classical orbits in the real space. 
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I. INTRODUCTION 



^ The fascinating properties of graphene originate from its 
"perfect two-dimensional character and from its unusual elec- 
C~-tronic structure where valence and conduction bands touch 
Tit two inequivalent points of reciprocal space. In the vicin- 
ity of these two points, the dispersion relation is linear 
■—a nd charge carriers behave as massless Dirac fermions 1 with 
striking consequences on the electronic properties, like the 
,dinusual 2D quantum Hall effect PS Recent studies have also 
addressed electronic properties of confined graphene struc- 
tures like dots, rings or nanoribbons. In particular, nanorib- 
3>ons have been suggested as potential candidates for re- 
_r>lacing electronic components in future nanoelectronic and 
djianospintromc devices! SH f n th a t perspective, the role of 
Sdge states is essential. In the absence of magnetic field, 
ledge states can emerge for particular types of edges.- 7 Such 
localized electronic states can be described analytically at 
Htow energy by using the Dirac equation with the appropriate 
7T)oundary conditions. 8 More recent theoretical works have in- 
i—aestigated the metallic nature of such edge states by consid- 
ering staggered on-site potentials, 9 surface perturbations p2 
T 7spin-orbit coupling,^ edge doppingP^l and sophisticated ter- 
J^iinations of the honeycomb lattice QMlSl Experimental ef- 
^Tprts have been performed to observe different natur es of 
O^dges^SHU] anc | study charge transport in ribbons!- 19 * 20 ! 
CN This paper deals with the structure of edge states in a 
r--»s.trong magnetic field, a problem which has been the sub- 
f^jfict of recent interest. Although this probl em is somehow 
(^Reminiscent of the case of usual 2D electrons] 21 * 22 ! it is much 
< — richer since here charge carriers are massless particules and 
^he structure of the edge states depends of the nature of the 
• 'Tjdges. In particular, Brey and Fertig first calculated within 
r%he tight-binding model the band structures of ribbons in a 
tjliagnetic field with two types of boundary conditions: the 
so-called zigzag and armchair edgesP^In addition, they have 
shown that these two types of edges can be described at low 
energy by a Schrodinger Hamiltonian with a harmonic po- 
tential. In the zigzag case, this potential is similar to the one 
of the usual integer quantum Hall problem, 21 whereas in the 
armchair case, the potential presents a specific asymmetry. 
Moreover, these edge states have been recently investigated 
numerically in graphene ringsj^and in graphene ribbons by 
considering an anisotropic hopping parameter,* 2 -^ or by taking 
account for a possible quantum Hall ferromagnetism! 2 ^ 
However, the energy dependence of the Landau levels 



when approching the edges has only been discussed briefly. 
The goal of this paper is to investigate in details the fine 
structure of the edge states energy spectrum for zigzag and 
armchair ribbons and provide a complete quantitative analyt- 
ical description within the WKB approximation. The paper 
is organized as follows: In section [Hi we recall and discuss 
in details the tight-binding spectra for zigzag and armchair 
ribbons with and without magnetic field. When a field is ap- 
plied, we point out peculiar features of the edge states which 
we describe analytically in the rest of the paper. Follow- 
ing the procedure suggested by Brey and FertigJ^S we derive 
in section III a simple effective Hamiltonian (essentially the 



squared Dirac Hamiltonian) with a potential depending on 
the boundary conditions. Then, in section |TV] we present a 
semiclassical framework to calculate analytically the low en- 
ergy spectra. Two methods are used: the first one is based 
on the Bohr-Sommerfeld quantization of the action and leads 
to a simple qualitative picture in terms of classical skipping 
orbits in section |Vj However, this method does not prop- 
erly describe the case where the classical cyclotron radius is 
of the order of the distance to the edge. Therefore, we de- 
velop a more sophisticated semiclassical approach, based on 
the WKB approximation, that leads to analytical quantita- 
tive results in agreement with numerical calculations. The 
details of this approach are given in the appendix. We con- 
clude in section IVII1 



II. TIGHT-BINDING SPECTRUM 

A. Zero magnetic field 

We briefly recall the band structure of infinite graphene 
ribbons within the tight-binding picture. For that purpose, 
we first consider the Hamiltonian H of an infinite 2D sheet 
of graphene, and we use the Bloch theorem in both x and y 
directions. As the honeycomb lattice has two carbon atoms 
per cell (that we call A and B, see Fig. IT]), the Bloch wave 
function is written as : 



1%) = 



V^ JkRj 
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(1) 



where Rj are the vectors of the triangular Bravais lattice. 
Then, within the tight-binding model, the Hamiltonian is 



written as : 



H = t 



f(k) 
f*(k) 



(2) 



in the basis of the two sublattices (^a,^b) where t is the 
hopping parameter between nearest neighbors and with : 



the two points K and K' for zigzag edge is not captured 
by the projection of the bulk result. This zero energy state 
is localized near the edges and its existence depends on the 
boundary conditions! 7 ! 8 ! 28 ! This edge state exists whitout any 
magnetic field, contrary to the edge states in the quantum 
Hall regime which are studied in this paper. 



f(k) = —1 — exp (ik ■ Si) — exp (ik ■ 0,2) 



(3) 



where a\ and 0,2 are the basis vectors of the triangular Bra- 
vais lattice with ||ai|| = \\a2W = clq. The lattice parameter 
a Q is related to the carbon-carbon distance a = 0.142nm by 

a Q — \/ia. The dispersion relation given by e(fc) = ±t f(k) 

(see Fig. pja)) consists in two bands which touch at the cor- 
ners of the first Brillouin zone (FBZ). The positions of these 
points are given by the condition f(K) = 0, that is : 



K = 3K -0*2) 



Z' = -\W 



where a\ and a* 2 are the basis vectors of the reciprocal lattice. 
The spectrum is linear near these points. The choices of the 
unit cell and of the axes are illustrated in Fig. fTVa) . The 
positions of the so-called K and K' Dirac points are shown 
in Fig. nib), where we have also indicated their projections 
onto the axis fc T and k„. 
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FIG. 1: (a) Sheet of graphene with zigzag and armchair edges, 
(b) The black dots represent the nodes of the triangular reciprocal 
lattice whose a* and a| are the basis vectors. The Dirac points K 
and K' are located at the corners of the FBZ. The projections of 
the K and K' points on the k x and k y axes have been represented 
by crosses. 

We consider now infinite ribbons with either zigzag or arm- 
chair edges. In these finite geometries, the Bloch theorem 
can be used only along the infinite direction, and the bound- 
ary condition along the finite direction yields a finite num- 
ber of bands. These bands, computed in the tight-binding 
model, are displayed in Figs. [2] (d-e). It is seen that the 
spectra of these ribbons correspond to the projection of the 
2D spectrum along the k x axis (armchair edge) or the k y axis 
(zigzag edge). The projected points K and K' coincide in 
the armchair case but they do not in the zigzag case. This 
remark illustrates the fact that for the armchair case, the 
two valleys are admixed by the boundary condition as we 
will discuss later. Finally, note that the flat level between 
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FIG. 2: (a) Energy spectrum e(k) of an infinite sheet of graphene. 
Projection of this spectrum (b) along the k y axis onto the plane 
(e, k x ) and (c) along the k x axis onto the plane (e,k y ). Band 
structure for an infinite armchair ribbon (d), and for a zigzag 
ribbon (e). Their widths are L = 49ao/2 for the armchair ribbon 
and L = 49y / 3ao/2 for the zigzag ribbon. The positions of the 
K and K' points are represented by dashed vertical lines and the 
ID FBZ is delimited by continuous vertical lines. The energies 
are given in units of the hopping parameter t. 



From now on, we define the geometry of the armchair and 
zigzag ribbons as shown in Fig. [3] : the ribbon is infinite 
along the y direction and it has a finite width L in the x 
direction. Within this convention, we have K y ^ K' for the 
zigzag ribbons and K y = K' for the armchair ones. 



B. Perpendicular magnetic field 

We now apply a perpendicular magnetic field B = —Be z 
on each type of ribbon and calculate numerically the band 
structure within the tight-binding model. 26 The finite width 
of the ribbons along the x axis leads us to work with the 
Landau gauge A y = —Bx. Because of the vector potential 
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FIG. 3: Graphene ribbons of width L with (a) zigzag edges and 
(b) armchair edges. We consider the y direction as infinite. The 
circles represent empty sites and define the edges. 




FIG. 4: Tight-binding spectrum at low energy for a zigzag ribbon 
with (red) and without (green) magnetic field. The dimensionless 
magnetic flux is <f> = 0.00126<^o an d the width is L — 199-^/300/2. 



A y , the hopping parameter t (taken as unity in the rest of 

the paper), takes a Aharon ov-Bohm phase t — > te 4 *o ■* ' , 
where 0o — h/e is the magnetic flux quantum. We now 
study the low energy spectrum for a given magnetic flux <fi 
through an elementary plaquette. In Fig. [4] we show the 
low energy tight-binding spectra for the zigzag ribbon with 
and without magnetic field in the same plot. The spectra 
for the armchair case are shown in Fig. [5] The zero field 
cones are progressively transformed into flat Landau levels 

with the expected e n = ±t\/2Tr\/3 n 4>/4>o behaviour P3 

We now comment on several important features of these 
spectra which to our knowledge have not been discussed in 
the literature. Since in the chosen Landau gauge, a state k y 
is centered at the position k y £ 2 B along the x direction, where 
£b = y/h/eB is the magnetic length, the variation of the 
energy levels can be interpreted as a function of this posi- 
tion, as redrawn in Figs. [6] and [7J In these figures, we have 
indicated the position of the edges along the x direction, at 
low energy. Indeed, the position of one edge is fixed by the 
position of the Dirac points in zero field and the other edge 



FIG. 5: Tight-binding spectrum at low energy for an armchair 
ribbon with (red) and without (green) magnetic field. The di- 
mensionless magnetic flux is cf> = 0.00126<?!>o and the width is 
L = 199a /2. 



B where L is the width 

(') 



is located at a distance Aq y = L 

of the ribbon and where q y — k y — Ky 1 . As the two val 
leys are not admixed in the zigzag case, this operation must 
be performed for both K and K' valleys, as seen in Fig. 
[6] Therefore, it is clear that the dispersive character of the 
levels corresponds to the edge states, as first discussed by 
Halperin in the context of the integer quantum Hall effect 
of massive particules in two dimensional gasesP3 However, 
the band structure for graphene is more complex and de- 
pends on the type of edge as we discuss now. For the zigzag 
edge, we first emphasize that, inside a given valley, the spec- 
trum is not symmetric: the energy levels are not identical on 
left and right sides. Moreover, exactly at the left and right 
edges, the energy levels ( edge take peculiar values : they al- 
ternatively take the value of higher bulk Landau levels e bulk . 
This correspondence is displayed by horizontal color lines in 
Fig. [6| More precisely, for the K valley, we have the rela- 
with n > (blue lines) on the left edge 



tion et 



-bulk 
fc 2n+l 



whereas we have €t 



-bulk 
~2n 



(green lines) with n > for 
the right edge. In addition, we have the same structure in 
the K' valley, provided the role of the two edges is permuted. 
This remarkable distribution of the edge states will be ex- 
plained in section |IVC| The armchair case, in Fig. [7] is 
quite different. Each Landau level is doubly degenerate, a 
direct consequence of the valley admixing. The degeneracy 
lifting takes place close to the edge with a non monotonous 
behaviour for one every two levels, that will be explained in 



section IV C In the following we present a qualitative and 
quantitative description of these spectra within a semiclas- 
sical analytical approach. 



III. LOW ENERGY EFFECTIVE HAMILTONIAN 

A. The bulk Hamiltonian 

At low energy, the Hamiltonian pi can be linearized 
around the two Dirac points K and K' . Expanding k as 
k = K^' + q in each valley, the tight-binding Hamiltonian is 




eA y = q y — eBx. 



FIG. 6: Tight-binding low energy levels in the quantum Hall 
regime as a function of the position k v £% along the x axis, for a 
zigzag ribbon. The edges are represented by vertical black lines. 
Green and blue horizontal lines indicate the position of the bulk 
Landau levels for comparison with the position of the levels at 
the edges. Note the mirror symmetry of the spectrum between 
the two valleys K and K' . The energy is given in units of t. 




FIG. 7: Tight-binding low energy levels in the quantum Hall 
regime as a function of the position k v £%, for an armchair ribbon. 
The edges are represented by vertical black lines. The energy is 
given in units of t. 



then replaced by a 4 x 4 linearized Dirac-like Hamiltonian : 



iq x - q y 





iq x -Qy 








iq x + q y 





-iq x + q y 



which describes the two uncoupled valleys in the basis 

{<PAi fBi ^'ai ^p'b)i where 7 = ^t and ao is the lattice spac- 
ing. The eigenfunctions , 3/ ± can be expressed as the super- 
position of the contributions of the two valleys, so that we 
can write : 



tf^r) 



JK 



Kr (<PA{r) 



±e 



<>? (<Pa(?) 
<f' B (?) 



iK> 



(5) 



In the magnetic field B = —Be zi and within the Landau 
gauge A y — —Bx, we perform the Peierls substitution q y — > 



Then we introduce the dimensionless 
variables x/£g ~> x, and x c — q y £s- Since q x = —id x , the 
Hamiltonian in a magnetic field reads : 



jao 
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d x - 


-(x- 


-X c ) 




(x- 


-X c ) 





J 





-d x 



Next, it is useful to work with the squared Hamiltonian H 2 
which is diagonal. We introduce the dimensionless effective 
Hamiltonian Heff as : 



H' 2 



1 1 

l B 



which defines four effective Schrodinger equations 



with 



u 



eff 



2 x 



Heff* = E n -S> 



V u (x) 

V d {x) 

V d (x) 

V u {x) 



(7) 



(8) 



(9) 



where the potentials V u (x) and Vd(x) are given by 



V u {x) = ^{x-x c ) 2 + - 
V d (x) = \{x-x c ) 2 -\ 



(10) 
(11) 



and are displayed in Fig [8] Thus, the effective Hamiltonian 
H e ff simply describes only two different harmonic oscilla- 
tors. Because of the non-commutativity of the operators x 
and d x , the two harmonic potentials V u (x) and Vd(x) are 
shifted by one energy level. From ^ and (pi, one obtains 
two different Schrodinger equations, one with the potential 
V u (x) for ip a and <p' B , and the other one with the potential 
Vd(x) for ip B and tp' A . 

The eigenvalues e„ of the original Hamiltonian ([7| are ob- 
tained from the dimensionless eigenvalues E n of the effective 
Hamiltonian H e ff as : 



7 a o 



*~B 



(12) 



In the Landau gauge we have used, the wave function is a 
plane wave along the y direction and now reads : 



*±(f) = 






(r) 



1 

B 

e iq v y 



Kt (<Pa{x)\ , ifr? (<p' A ( x ) 



Vb{x) 



±e l 



ip' B {x) 



(13) 



The eigenvalues E n = n are the Landau levels associated 
with the eigenfunctions of the harmonic oscillator {$„} as 
follows : 



ifiA = $n-l <Pa = $n 
<PB = ®n <f' B = -$n-l 



(14) 




FIG. 8: The harmonic potentials V u and Vd of the effective Hamil- 
tonian Heff- The eigenfunctions associated with the eigenvalues 
E n = n are the eigenfunctions of the harmonic oscillator, respec- 
tively $„_i and <!>„. 



This is illustrated in Fig. [8] The effect of the edges is to mod- 
ify the potentials ( 10 1 and (Till) and consequently the compo- 



S>) 



nents f A / B of the wave functions and the energy spectrum. 
These modifications depend on the nature of the edge which 
yields to specific boundary conditions. In the two following 
subsections, we derive the effective Schrodinger equations for 
zigzag and armchair edges. 



B. Zigzag edges 

To treat the zigzag edges of a graphene ribbon, we first 
recall that all the atoms on one edge belong to the same 
sublattice, and therefore, all the atoms on the opposite edge 
necessarily belong to the other sublattice. As the ribbon is 
periodic in the y dire ctio n, we still write the low energy wave 
functions i S ± as in (13k but the components <fiJ/ B are not 
the eigenfunctions oi the harmonic oscillator anymore. As 
the left edge is only made of A sites (see Fig. |3|, the wave 



function on the B sites vanishes for x = i.e. ^^{x = 0, y) = 
0. The situation is identical on the right edge at x = L for 
the A sublattice : ^ A (x = L,y) — 0. As a consequence, we 
have the boundary conditions : 



<p A (L) = (V A (Z) = 
<p B (0)=0<p' B (0) = 0. 



(15) 



These boundary conditions do not admix the valleys which 
can still be described separately, as in the bulk system. Con- 
sequently, with the zigzag boundary conditions, we have four 
independent Schrodinger equations, one for each valley and 
for each sublattice, with four independent potentials Va(x), 
Vb(x), V a (x) and V B (x). Such boundary conditions can be 
accounted for by an infinite potential barrier either at x = 
or at x = L. acccording to the sublattice and the valley. 
Namely, the potentials V B (x) and V B (x) are harmonic po- 
tentials that have to be cut in x = whereas the potentials 
Va{x) and V A {x) have to be cut in x = L. Since the four 
edge problems are quite similar, we focus on the case of the 



right edge for the K valley. Thus, the potential Va(x) reads : 
V A (x) = { K(X) *» X<L (16) 



and is plotted in Fig. M Similarly, one can define V A (x), 
V B (x) and V B (x). As a matter of fact, the problem of a 
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FIG. 9: Potential Va(x) that describes the right edge in the K 
valley of a zigzag nanoribbon for the 1-Le.ff problem. The edge is 
modelized by an infinite potentia barrier. 

harmonic well with an infinite potential barrier (Fig. [9]) is 
identical to the problem of a double symmetric harmonic 
well (Fig. 10), provided we only keep, in the latter, the 



eigenenergies associated with the eigenstates that vanish on 
the edge, that are the antisymmetric states. The interest of 
considering such a double symmetric well is that, as we will 



see in section IIIC both zigzag and armchair edges can be 
described by double harmonic wells, the difference being an 
energy shift between the two wells in the armchair case. We 
illustrate in Fig. 10 the double symmetric potential Va(x) 



we finally consider in the effective Schrodinger equation for 
the right zigzag edge and the K valley, that is : 



1 



d 2 x + V A {x) >$>as(x) = E^(x c )^ AS (x) (17) 



V A (x)= l -{\x\+x c f + \ 



(18) 



The index AS refers to the antisymmetric solutions and we 
have performed the translations x — L — > x and x c — L — > x c , 



so that the right edge is now located in x = 0. The eigenen- 
ergies E^ s (x c ) are calculated analytically within a WKB 
treatment introduced in section lTVl The detailed calculations 
are explained in the appendix. 



functions now read 




X\ x c X 2 



FIG. 10: Double symmetric harmonic potential with x c = L — 2. 
We show the eigenenergies obtained within the WKB approxi- 
mation. The higher energy levels (blue continuous lines) are the 
antisymmetric solutions. They are the same than those of a har- 
monic well cut by an infinite potential barrier (Fig. [91. For 
completion we have indicated the energies of the symmetric solu- 
tions (dashed red lines) which are nore considered here. The two 
regions A and B and the cyclotron radius R are defined in section 
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C. Armchair edges 

We now consider the case of an armchair ribbon as dis- 
played in Fig. [3jb) . The corresponding low energy Hamil- 
tonian is obtained from (H| by the substitution q y — > —q x 
and q x — > q y . By keeping the same gauge A y = —Bx as 
above, the Hamiltonian in the presence of a magnetic field 
now reads : 
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-d x 
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(19) 



Of course, this new choice of axes does not affect the eigenval- 
ues of the bulk problem and leaves the effective Hamiltonian 
H e f f defined in (|7| unchanged, but the eigenfunctions have 
now the followingform : 



tp A = -i$ n _i (p' A = i$„ 

tpB = *« <Pb = ~$n-l 



(20) 



where the {$„} functions are the eigenfunctions of the har- 
monic oscillator. The phase factor i originates from the ori- 
entation of the ribbon. Following the procedure suggested 
by Brey and FertigJ^ we now construct an effective Hamil- 
tonian which accounts properly for the armchair boundary 
conditions. Since, for this type of edge, we have the par- 
ticular relation K y — K , as argued in section III the wave 



*±(r) 






J(K y +q y )y 



(21) 



AK x x 



,w± 



ip B (x) 



±e 



iK'x 



(f'Ux) 
ip' B (x) 



where the components <f B / A must satisfy specific bound- 
ary conditions. A crucial difference with zigzag ribbons is 
that armchair edges are constituted of both A and B sites. 
As a consequence, both components ^ A and \& B must van- 
ish on each side of the ribbon: i \! ± (x — Q,y) — and 
\l/ ± (a; = L,y) = 0. From now on, we choose to explicit 
these boundary conditions for the left side at x — 0, what 
leads to : 



vi(o) = 



T^(0). 



(22) 



Contrary to the zigzag case, the armchair boundary condi- 
tions admix the contributions of the two valleys for a given 
sublattice. In addition, since the Dirac Hamiltonian is first 
order, the continuity equations (22 1 on the wave function 
imply continuity of the derivatives, 
obtain : 



From (191 and (22 1 we 



d x Va 
d x (p B 



±d x v'£ 

±d xV 'i 



(23) 



Next, as suggested by Brey and Fertigp^ we build new func- 
tions <& ± as : 



$±(x) 



tti(-x)e(-x)T<f£(x)6(x) 

-i[cp't(-x)e(-x)T<fi(x)e{x)] 



(24) 



where 9 is the Heaviside function. Thus, the functions <i> ± 
are solutions of a new effective Schrodinger equation with 
a potential Vi e ft(x) which is V u (x) for x > and Vd(x) for 
X < : 



1 



d x + V left (x) )<I> ± (x)=E±(x c )<S> ± (x) 



Vleft{x) 



\x\ - x c ) 2 + 6{x) - 6{-x) 



(25) 
(26) 



In ord er to t ake into account the armchair boundary condi- 



tions ( 22JI23 1 , we impose the functions $^ /B to satisfy the 



A/B 



continuity equations : 



$ ± (o+) = $ ± (o-) 



dx<S>~- 



0+ 



= fl«* ± 



(27) 



The asymmetric potential (26 1 is shown in Fig. 11 The 
asymmetry originates from the non-commutativity of x with 
d x and corresponds to a shift in energy by one Landau level 
between the two uncoupled wells V u and Vd- The asymmetric 
structure of Viejt is a consequence of the valley admixing 
imposed by the boundary conditions. 

We now discuss the qualitative structure of the wave func- 
tions. The two combinations ^> ± are built with four compo- 
nents (p A / B - When the particule is far from the edge, the 
system is described by two independent potentials V u (x) and 
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FIG. 11: Potential V; e /t(a;) for the left edge of an armchair ribbon 
in a magnetic field. The two wells are centered from a distance 
\x c \ to the edge x — 0. The energy levels E% in blue and E~ 
in red are shown by horizontal lines in each well and have been 
calculated semiclassically. The energies are dimensionless. 



Vd(x) so that the components <pJ/ B take their bulk value 
given by (20 1. Now, close to the edge, the expressions of 
In particular, the boundary conditions 

by matching 



fjl/B are mo dified 

(27 1 imply two types of solutions $ + and $ 



the components <pJi B in two different ways. We illustrate 
this point in Fig. [l2]in the case where n — 1 by showing a 
qualitative construction of the components (fi A / B - The eigen- 
function $ + has In nodes whereas <I>~ has 2n — 1 nodes. As 
a consequence, the associated eigenvalues E~(x c ) are lower 
than the eigenenergies E+(x c ), and the valley degeneracy of 
the energy spectrum is lifted, as seen in Fig. [5] This point 



is discussed quantitatively in section IV C 

By imposing ^f(L, y) = 0, the right edge can now be con- 
sidered easily in the same way. We obtain a similar effective 
problem where the potential Vi e ft(x) has been replaced by : 



V„ght(x) 



1 



(M 



■ e(-x)-e(x) 



(28) 



FIG. 12: Illustration, for n = 1, of the two solutions &~ (x) and 
$ + (a:) at the left armchair edge. These solutions do not exhibit 
the same number of nodes (black circles), (a) The continuity 
equations generate a wave function <I> _ with 2n — 1 nodes. The 
fonctions tpg (—x) and —ifp'^ {—x) are represented on the left (red) 
whereas ip'g (x) and —itp~^(x) are represented on the right (blue), 
(b) The continuity equations generate a wave function $ + with 
2n nodes. The functions <£%{— x) and — iip'^(— x) are represented 



on the left (red) whereas (p'g (x) and 
the right (blue) 



iip\(x) are represented on 



\j' rightyX) 




4 


E 






^V 




3 


i 


\ 


I 


\ / 


: / 


T> \ 




/• 






/ 


\, Ri 


. 


/ . 




, Rr 


/ 


C ^ i 




-4 Xl 


x c 


x 2 


i\ 

edge 


2 


t x 



where we have performed the translations x — L — > x and 
x c — L — > x c for more commodity in the calculations in the 
next section. Thus, the potential V r i g ht(x) that describes the 
right edge at x = L is now centered in x = 0. The potential 
Vright(x) is displayed in Fig. [l3 



FIG. 13: Potential V r ight(x) of the right armchair edge in a mag- 
netic field with the dimensionless energies E n calculated semiclas- 
sicaly. Here x c — —2. The three regions C, T> and £ , such as the 
cyclotron radii Ri and R r refer to the section [TV] 



IV. SEMICLASSICAL TREATMENT 

The aim of this section is to calculate analytically the en- 
ergy spectrum E n (x c ) for both zigzag and armchair bound- 
ary conditions, by using a semiclassical formalism. The 
eigenenergies are solutions of effective Schrodinger equations 
in appropriate potentials: a double symmetric harmonic po- 
tential Va{x) (Fig. 10) for the zigzag case and a double 



asymmetric harmonic potential V r i g ht(x) (Fig. 
armchair case. We develop here two approaches. 



13] for the 
In the first 
one (Sec. IV A), we assume that the two wells of the poten- 



tials are uncoupled, and quantize the action with the Bohr- 
Sommerfeld rule. The second approach (Sec. IV B), based 



on the WKB formalism, properly accounts for the overlap of 
the wave function between the two wells. 

From now on, we only focus on the right edge of the rib- 
bons, and take H = 1. 



A. Semiclassical quantization of the action 

We first define the classical cyclotron radius R. For the 

zigzag case, R is given by E = l/2(i? 2 + 1) as illustrated in 

For the armchair case, since the potential V r i g ht(x) 



Fig. 10 



is asymmetric, we need to define two classical cyclotron radii 
Ri and R r (see Fig. 13 1. Next, we introduce an action S 
associated with these potentials. Along a closed path, the 
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action is given by : 

S(E,x c ) = 2 / dx y/2{E-V(x)) (29) 



where xi and xi are the positions of the turning points, and 
V(x) is either Va(%) or V r i g ht{x). Such an action is quantized 
with the Bohr-Sommcrfcld rule as : 



S(E,x c ) = 2ir(n + j(E,x c )) 



(30) 



where n > is an integer and < j(E, x c ) < 1 is a function 
that encodes all the information on the connection procedure 
at the turning points. For a harmonic potential. 7 = 1/2. 
Here, for the double well potentials, 7(x c ,_E) is not a con- 
stant anymore, but depends on the energy and on the dis- 
tance to the edge because of the overlap of the wave function 
between the two wells. We assume in this section that the 
two wells are uncoupled, so that we take 7 as a constant. 

Now we have to specify distinct regions in energy, each 
of them requiring an appropriate semiclassical treatment. 
These regions are delimited by horizontal thick lines in Figs. 
[TO] and [13] and involve different expressions of the action 
S(E,x c ) given by (29} . 



1. Zigzag ribbons 

We consider the potential V A (x) that describes the right 
zigzag edge for the K valley (Fig. 10 1. We recall that the 
right edge is located at X = 0. For a given energy E, we 
have to discuss two distinct regions : \x c \ > R (region A) 
and \x c \ < R (region B), where R = \J2E — 1 is the cyclotron 
radius. 

• Region A: \x c \ > R. The turning points are defined by 
X\ = x c — R and X2 = x c + R so that the action (29 1 
in the left well is simply : 



S A (R) = irR 2 



(31) 



A simple calculation of the energies by quantizing S A 
with the constant value 7 = 1/2 in (30 1 leads to the 
Landau levels E n = n + 1 with n > 0, and therefore to 
the spectrum e n ex y/B(n + 1) in this region. The n+1 
term originates from the energy shift of the potential 
Va(x). The rest of the spectrum near the right edge is 
given by the contribution of the other valley. There- 
fore, by treating the same way the potential V A (x), 
we find the spectrum E n = n with n > 0. Finally, 
we obtain a set of degenerated energies E n — n with 
n > 1, and a non degenerated level Eq = what gives 
the expected valley degenerated graphene energy levels 
e n ex zkyBn for n > 0. 

Region B: \x c \ < R. The left turning point X\ is 
unchanged but the right turning point is replaced by 
X2 = so that the action reads : 



.S'tf(/?..r ( ) = If (^-^sin(20) 



(32) 



where we have introduced the parameter = arccos %■ . 
The spectrum depends now on the distance x c to the 
edge. The total action in the double well St = 25$ is 



quantized as St = 2n(p + 1/2). We recall that we have 
to keep only the antisymmetric solutions which have 
an odd number of nodes, what implies p = 2n + 1. 
Therefore, such a quantization leads to : 



S B {R,x c ) = 2ir 



(33) 



We recover the result 7 = 3/4 for a harmonic potential 
cut by an infinite barrier potential. 34 By identifying 
the two expressions ( 32 1 and ((33]) , we finally obtain a 



set of self-consistent equations labeled by the integer 
n for E n (x c ), from which we can extract the energy 
spectrum. 

The spectrum E n (x c ) obtained within this approach for 



the right zigzag edge in the K valley is plotted in Fig. 15 
In order to correctly describe the region where \x c \ ~ R, we 
use in section |IVB| a more sophisticated approach based on 
the WKB formalism. 



2. Armchair ribbons 

Because of the asymmetry of the potential V r i g ht(x), we 
have now to distinguish three regions (C, T>, £ in Fig. 13 1. 



In each one, we consider two actions, Si for the left well and 
S r for the right well. In addition, the levels in the left or in 
the right well will be indexed by different integers ni,n r > 0. 
The energy is given by E = (Rf + l)/2 = (Rf. — l)/2 where 
Ri (resp. R r ) is the cyclotron radius for the left (resp. right) 
well. 



In each well, the action still 



Region C: Ri < R r < 

has the form given by (31 1, so that we write : 



left well : Si = S A (Ri) 
right well : S r = S A (R r ) 



(34) 



By using the Bohr-Sommerfeld rule (30 1 with 7 = 1/2 
for both Si and S r , we find the energies E ni = m + 1 
into the left well and E nr — n r into the right well, 
and then n = n r = n\ + 1. This leads to the degen- 
erated Landau levels E n = n with n > 1 and a non 
degenerated level Eq = 0. 

• Region T>: Ri < \x c \ < R r . Now, because of the step 
in the potential, the action S r has the form given by 
( 32 1 , whereas Si is unchanged. The two actions read 



left well : Si = S A (Ri) 
right well : S r = St3(R ri x c ) 



(35) 



The action Si is still quantized as previously, what sim- 
ply gives the spectrum E n — n with n > 1. As the wave 
function inside the right well does not totally vanish in 
X = 0, the situation is different than in the case of an 
infinite potential barrier. This suggests a choice closer 
to 7 = 1/2 rather than 7 = 3/4 for the quantization 
of the action S r . By doing so, we obtain an implicit 
equation in E n and x c from which we can extract the 
spectrum E n (x c ). Of course, we understand that the 
region T> is necessarly badly described by such an ap- 
proach and specifically requires a more sophisticated 
method since 7 cannot be a constant anymore. 



• Region 8: \x c \ < Ri < R r - In this region, we need to 
consider the total action that reads : 



Ss = Sb(Ri,x c ) + S B {R r ,x c ) 



(36) 



where the function Sb is given in (p2|. This expres- 
sion is still valid when x c crosses the edge. Since the 
two turning points both touch a harmonic potential, 
this action is simply quantized by the usual Bohr- 
Sommerfeld rule with 7 = 1/2 what also leads to an 
implicit equation in E n i and x c , where n! is a different 
integer than n, since there are twice more solutions in 
this region than in the regions A or C. 

The armchair spectrum obtained within this approach is 
shown in Fig. 
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It perfectly matches with numerical so- 
lutions of the Schrodinger equation ( |25| except around the 
intermediate region T>. In the following, we describe quanti- 
tatively all the regions by keeping 7 = j(E, x c ) and calculat- 
ing the spectra of zigzag and armchair ribbons by accounting 
for the coupling between the two wells. We give now the gen- 
eral picture of the method which is detailed in the appendix. 



B. The WKB approach 



the potential V(0). The same procedure must be performed 
for the right well where we obtain similar relations between 
the ratio a r //3 r and the action S r calculated between the 
two turning points in the right well : 



13 ^ : tanSV/2 = a r //3 r 

X3 = : tan(SV/2 + S r ) = a r //3 r 

where x% is the left turning point in the right well. 



(41) 




Xl X c X2 

(a) 




(b) 



In this section, we present a method based on the WKB 
formalism to account for the overlap of the wave function 
between the two wells. This approach is detailed in the ap- 
pendix. Technically, it brings us to calculate the function 
"f(E, x c ). For this purpose, we express the wave function 
inside the left well within the WKB approximation as : 



^ {x) = 7fcp sin ( 5(:El,;E) + i) 



with C a constant and where 



S(xi,x)= / dx'y/E - V(x') 



(37) 



(38) 



is the partial action between the turning point x\ and an 
arbitrary position x inside the left well. Close to the right 
turning point X2, this approximation breaks down. In or- 
der to find a valid approximation of the wave function in 
this region, we linearise the harmonic potential so that the 
cigcnfunctions of the Schrodinger equation are a combination 
of Airy functions, that is : 

<p-+(x) = a,Ai(/,(a:, E, O) + ftBi(/,(s, E, x c )) (39) 



where ai and /3i are constants and fi(x, E,x c ) is a function 
that depends on the region A, B, C, T> or £, so that each 
region must be treated separately. By imposing the equal- 
ity between ip^(x) and the asymptotic expansion of ip~*(x) 
inside the well, we obtain a relation between the constants 
and the action Si which is given by the relation ( [29| . We 
need to distinguish the cases X2 ^ and X2 = (Fig 
for which the matching condition gives : 



14) 



x 2 ^0: tenSi/2 = cn/pi 

x 2 =0: tim(Si/2 + 8i)=an/Pi 



(40) 



FIG. 14: Illustration of the two typical different situations we 
have to distinguish for the semiclassical calculation, (a) The two 
turning points x\ and X2 are on the parabola, (b) The turning 
point X2 = is not on the parabola. These two situations involve 
different expressions of the action S(E, x c ) = Si in the left well. 
Same considerations have to be made for the right well. 



To determine the ratio ai/ Pi, we impose the continuity of 
the wave function and its derivative between the two wells. 
By doing so, we obtain one equation with the two unknown 
quantites ai//3i and a r / fi r . Next, by injecting S r = St — Si 
into (41 ), where the total action St is calculated explicitly as 
we obtain a second order polynomial 



a function of E and x c 
in ai/Pi what yields two solutions. Because of the relation 
(40), the action is quantized as Si — 2n(n + 7/) which is the 
Bohr-Sommerfeld quantization. Thus, the two solutions of 
the polynomial yield the two solutions Ji(E, x c ). 



where Si = Si(E, x c ) is a known function originating from the 
step between the energy E at the turning point X2 = and 



Next, by using the formula (29 1, we calculate explicitly the 
expression of the action Si as a function of E and x c , which 
we identify with the action given by the Bohr-Sommerfeld 
rule previously found. As we know the functions r ) l (E,x c ), 
we finally obtain two analytical self-consistent equations in E 
and x c for a given n, one with 7^ , the other one with jf . The 
spectrum E n (x c ) can then be extracted for a given n. The 
whole procedure must be performed for each region A to £ 
where the action has different expressions as discussed in the 
previous section. All this study is detailed in the appendix 
for both the symmetric and the asymmetric potentials. 



C. Quantitative analytical results 

We give here the spectra obtained for the effective poten- 
tials Va{x) (zigzag case (18)) and V r i g ht(%) (armchair case 
(28)), within the two semiclassical approaches introduced 
above. 
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1. Symmetric potential (zigzag) 

The energy levels E n (x c ) of the double symmetric har- 
monic potential Va (x) describing the contribution of the K 
valley at the right zigzag edge (see Fig. 10 ) are shown in Fig. 
15 The two regions \x c \ > R and 
a parabola of equation E — (x 2 4 
Fig. 15 we have only kept the higher energy solution of the 



x c \ < R are separeted by 
l)/2. In the spectrum of 



double well problem, since they correspond to the eigenener- 
gies of the antisymmetric eigenfunctions. The dashed lines 
correspond to the approximation where 7 = est, what leads 
to an unphysical discontinuity around \x c \ — R. The con- 
tinuous lines represent the energy levels obtained within the 
WKB approximation. The solutions found with this method 
in regions A and B perfectly match at \x c \ — R. Next, we 



\ B / / 

\ s / s 

A \^S / 






E(x c )/ 








right edge 










FIG. 15: Dimensionless eigenenergies E n (x c ) of the double sym- 
metric harmonic oscillator potential Va(x) calculated analytically 
(continuous line) within the WKB method and (dashed line) 
within the Bohr-Sommerfeld quantization rule with 7 = | in re- 
gion A and 7 = I in region B. The solutions we keep correspond 
to the antisymmetric solutions of the double symmetric well. 



show in Fig. [16] the structure of the edge states spectrum 
in the K valley for both the left and right edges. The right 
(blue) part of the spectrum, is the energy spectrum on the 
right edge calculated within the WKB method and already 
displayed in Fig. 



1G 



whereas the left (red) part represents 
the edge states on the left side of the ribbon in the same K 
valley. The effective problem for the left edge is obtained in a 
similar way than for the right edge (see section III B ) where 
we have considered the potential Vd{x) and the boundary 
condition (Pb(x = 0) = 0. Note that we recover the partic- 
ular relation eff 9e = e^n+i with n > on the left edge and 
£ edge _ £ buik w j^ n 77, > for the right edge. This peculiar 
property of the edge states is due to the fact that we have 
to keep only the antisymmetric eigenfunctions of the dou- 
ble well problem. In particular, when x c — 0, the double 
well is simply a single harmonic well, so as the bulk pote- 
nials. In addition, the shift of this distribution between the 
two edges is due to the shift in energy between the poten- 
tials Va(x) and Vb(x) that respectively originate from the 
potentials V u (x) and Vd(x) modified by the zigzag bound- 
ary conditions. We have emphasized this point by dashed 
lines in Fig. WE\ This explains the remarkable structure of 



the zigzag edge states discussed in section [II B| Finally, the 
mirror symmetry between the valleys K and K' emphasized 
in Fig. |6]is now clear, since the edge potentials into the K' 
valley, that are V A (x) for the right edge and V B (x) for the 
left edge, now originate from the bulk potentials Vd(x) and 
V u (x) respectively with the zigzag boundary conditions. 




FIG. 16: Dimensionless zigzag edge states E n (x c ) of the Heff 
problem in the K valley. We have represented the bulk wave 
function of each sublattice (left (red) for B sublattice and right 
(blue) for A sublattice). The results have been obtained within 
the WKB approximation. 



2. Asymmetric potential (armchair) 

For the armchair case, we have to keep all the solutions 
of the double asymmetric potentials Vi e ft(x) and V r i g ht{x). 
As discussed in section [III C| and illustrated in Fig. [12] the 
wave functions of the effective Schrodinger equation with 
armchair boundary conditions ( 25 ) have, for a given n, either 
2n nodes or 2n—l nodes, what leads to pairs of eigenenergies 
En(x c ). The semiclassical energies E n (x c ) are displayed in 
Fig. [IT] The three regions C, T> and 8 are separeted by 
the parabola parametrized by |x c | — Ri and \x r \ = R r , that 



are E — \{x 2 c + 1) and E = \(x* c — 1). Contrary to the 
simple semiclassical calculation with 7 = est (discontinued 
brown lines), the WKB method gives a perfect matching 
between the three regions and is in very good agreement 
with numerical calculations (grey crosses) . 



V. QUANTIZED SKIPPING ORBITS 

In this section, we propose a simple interpretation of the 
results in terms of a semiclassical picture for skipping or- 
bits. Consider the action S associated with the effective 
Schrodinger equations with the potentials obtained for the 
zigzag case. These four equations (one per edge and per 
valley) describe the motion of a free massive particle in a 
magnetic field (all units being set to one) in the presence of 
an infinite potential barrier or the image of the harmonic po- 
tential, as illustrated in Figs. [9] and [TO] for the right edge in 
the K valley. From the classical equations of motion, the ac- 
tion S = <f p- dr along a closed trajectory for a given energy 
is simply related to the area A enclosed by the corresponding 
periodic orbit : 



S(E,x c 



p ■ dr 



eB I - 1- 

— (b r x dr 

2 J 



eBA 



(42) 
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FIG. 17: Dimensionless energies E n (x c ) of the armchair edge 
states. The results obtained with three different methods are 
shown. The discontinuous brown lines are the levels obtained by 
taking 7 = est in the Bohr-Sommerfeld quantization (BSQ). The 
colored continuous lines have been obtained within the WKB ap- 
proach. The different colors in regions C, T> and £ are associated 
to the different semiclassical behaviours of the particule in Fig. 
|19| as discussed in the next section. The grey crosses represent 
the numerical results obtained within the tight-binding method. 
We show that the results obtained within the WKB approxima- 
tion perfectly fit the numerical calculations, except for n = 0. 



so that the Bohr-Sommerfeld quantization rule ( 30 ) implies 
the quantization of the area 



A(E,x c ) = 2-n(n- 



•7)4 



(43) 



in our units where h = 1. In the bulk of the system, that is 
when the distance |a; c | to the edge is larger than the cyclotron 
radius R, we have simply 



A(E) = S(E) = 2n{n + 7) 



(44) 



and the distances are measured in units of £b- This quantiza- 
tion implies the quantization of the cyclotron radius R. Tak- 
ing 7 = 1/2, A = 71- i? 2 = 2n(n + 1/2), so that R\ = 2n+l, 
we deduce the energy levels, E n = R^/2 ± 1/2 = n, and 
obtain the bulk spectrum of Landau levels E n = n. 

When \x c \ becomes smaller than R, the classical orbits 
skip along the wall represented by the infinite potential well, 
and they form open orbits. As already noticed 35 , the area 
to be quantized is the area delimited by the skipping orbit 
and the wall. In this case, the action calculated in (32 1 has 



a very simple interpretation in terms of skipping orbits since 
the parameter = arccos ^ has a clear geometrical meaning 
shown in Fig. [THJ The action is nothing but the area enclosed 
by the skipping orbit : 



A{E,x c ) = R 2 



-\sm(26) 



(45) 



We can now interpret semiclassically the energy spectra ob- 
tained in the above subsections. 



For the zigzag edge, the two regions A and B in the spec- 
trum (Fig. 15) can be associated to the two distinct cases 



represented in Fig. 18 

(a) region A, A(E) = irR 2 = 2n 

(b) region B 



(46) 



R 2 / 3 

A(E, x c ) = — (29- sin 26>) = 2tt I n + - 



from, which together with the relation E — R 2 j2 ± 1/2, we 
have deduced the energy levels, except near the region x c ~ 
—R. For the armchair case, the situation is more involved. 






edge 



(a) 



edge edge 

(b) 



FIG. 18: Classical cyclotron orbits for the "effective" particle 
associated to the effective Hamiltonian Heff close to a zigzag 
edge, (a) The closed cyclotron orbit encloses an area A quantized 
by the Bohr-Sommerfeld rule |47| a nd the resulting Landau levels 
spectrum is displayed in Fig. |15| (region .4). (b) The particle 
is skipping along the edge and the edge states energy spectrum 
(region B in Fig. 151 is obtained from the quantization rules (47 1. 



since we describe the classical motion of a massive particle 
in the presence of its shifted image (Fig. 19 1. Since the 



potential is asymmetric, we have defined two cyclotron radii 
related to the energy E as E = Rf/2 + 1/2 = R 2 /2 - 1/2. 
Now we have to consider three different cases corresponding 
to the regions C, T> and £ in Fig. 17 

In the region C where \x c \ > R r > Ri, the orbit and its 
image are both closed cyclotron orbits. The quantization of 
their area : 



left side: Ai = ttR? = 2tt [ n t + - 



right side: A r = kR t = 2tt [ n 



(47) 



leads to the the energy levels given by E n — n, with n = n r = 
m + 1, what quantize the cyclotron radii as R r = \j2n + 1 
and R r — \j2n — 1, the latest does not allow the n = level. 
In the region T> where R r > \x c \ > R[, the left orbit is 
still closed but its image is a skipping orbit. Keeping the 
notation n = n r = ni + 1, we have : 



left side: Ai = nRf = 2n(n - |) 



n >0 



right side: A r = R 2 



isin26>) = 2?r(n+ \) n > 



(48) 



from which we obtain the energy levels in region T>. The 
degeneracy lift is classically explained by the fact that both 
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orbits have a different structure. The left orbit is closed and 
its semiclassical energy is still E n = n + 1, while the energy 
associated with the right skipping trajectory increases with 
x c . Note that it is in this peculiar region that an inversion 
of the slope occurs. 

In the region 8 where R r > Bi > \x c \, both orbits are 
open and the area to be quantized is the total area 



Ai + A r = 2?r n' 



ri > 



(49) 



Since the action has almost doubled from region T> to region 
£ , the resulting energy spectrum is twice denser, and the 
two-fold degeneracy of the levels has been removed. 




edge 




X c -X c X 

(a) 





edge 
X c - X c 

(c) 




VI. TWO EDGES 

Up to now, we have considered the evolution of the spec- 
trum in the vicinity of one edge. If the two edges are suf- 
ficiently far apart compared to the cyclotron length lg, the 
spectrum can be treated independently on both sides. We 
now consider the case of a narrow ribbon, whose width L is 
of the order of a few magnetic lengths. The low energy spec- 
trum in a ribbon of width L — 9.6£b is shown on Fig. |20[a) 



and 21 b), respectively for zigzag and armchair ribbons. 




FIG. 19: Classical cyclotron orbits for the "effective" particle 
associated to the effective Hamiltonian, close to an armchair edge, 
(a) Each cyclotron orbit encloses an area Ai or A r , each of them 
quantized by the Bohr-Sommerfeld rule (30 1. For an energy level 
E n — n, one has Ri = \/2n — 1 and R r — yJ2n + 1. The ground 
state n — has no component on the left side. The associated 
spectrum is the bulk Landau levels spectrum displayed in the 
region C in Fig. |17[ b) When the left orbit approaches the edge, 
its image becomes a skipping orbit. This lifts the degeneracy 
between the energy levels as seen in Fig. [17] (region T>). (c) and 
(d) When x c is closer to the edge, both orbits are open skipping 
orbits. The quantization of the total area leads to the spectrum 
in region £ of Fig. |17| 



FIG. 20: (a) Red curves: Low energy spectrum for a zigzag rib- 
bon in a magnetic field in units of t. The width L of the ribbon is 
L — 174a and the dimensionless magnetic flux / is / = 0.00126, 
so that L/Ib = 9.6. The black curves indicate the situations 
where the classical cyclotron orbits touch an edge, (b) Schematic 
representation of the different situations for the position of the 
cyclotron orbits with respect to the edges of the ribbon, as dis- 
cussed in the text. 

In the zigzag case, in each valley, the spectrum clearly 
exhibits three different regions, which correspond to the ge- 
ometry of the orbits depicted in Fig. 
corresponds to bulk cyclotron orbits. 



20 » 



Regions A 
Regions Bi and B r 
have been already discussed and correspond to a skipping 
orbit along a single boundary (left or right). Each sublattice 
is characterized by a cyclotron orbit (in K valley, small blue 
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orbit for A sites and large pink orbit for B sites; the oppo- 
site for the K' valley). As explained in the text, in K valley, 
the A cyclotron orbits sees only the right edge and the B 
cyclotron orbits sees only the left edge (the opposite in the 
K' valley) . At high energy, the new region J corresponds to 
the situation where both cyclotron orbits intersect the two 
boundaries. The black curves correspond to the situations 
where a cyclotron orbit precisely touches a boundary. Their 
equation is 



,^0 



e(k y ) = t?-p-y/( qy -A qy )2±l/P B 



(50) 



where Aq y — for the left edge and Aq y = L/£ B for the 
right edge. By a simple quantization of the areas shown in 
Fig. 20 b), one can obtain the full low energy spectrum in 



a very good approximation, excepted when the orbits graze 
the boundaries, that is in the vicinity of the black curves. 
It is obvious on this figure that the two edges do not play 
exactly the same role, and that, for a given valley, the spec- 
trum is not exactly symmetric. Furthermore, in each valley, 
an additional dissymetry is seen at high energy because the 
linear approximation of the low energy Hamiltonian breaks 
7| down. 

The case of armchair boundary conditions is more in- 
volved. In Fig. 21 a), one can distinguish 9 regions, which 



(7) 



corresponds to the geometries of the classical orbits depicted 
in Fig. 21 'b). The black curves correspond to the situations 
where a cyclotron orbit precisely touches a boundary. Their 
equation is given by Eq. ( 50 I . Region C corresponds to bulk 



cyclotron orbits. The blue cyclotron orbits are related to the 
A and B' eigenfunctions, whereas the pink cyclotron orbits 
are related to the A' and B eigenfunctions. Regions £i and 
£ r have been already discussed and correspond to skipping 
orbits along a single boundary (left or right). In the interme- 
diate regions T)\ and T> r , only one of the two orbits touches 
the edge. At high energy, in the region I, the two cyclotron 
orbits touch the two edges of the ribbon. The more exotic 
regions J- . Qi and Q r have a simple geometric interpretation 
in Fig. |2lTb). 



VII. CONCLUSION 




We have investigated the spectra of graphene ribbons in 
a magnetic field with zigzag and armchair boundary condi- 
tions. We have first revisited these spectra numerically and 
revealed a remarkable structure in the repartition of the edge 
states in zigzag ribbons, that we explain in this paper. We 
notice that these remarkable behaviours of the edge states 
could be observed by scanning tunneling microscopy (STM) 
or spectroscopy (STS) techniques) 30 " 33 ! 

Next we have described and calculated these edge states 
with simple analytic tools. For both types of ribbons, 
effective Schrodinger equations with a specific double well 
potential have been derived at low energy. This potential 
is asymmetric in energy when the valleys K and K' are 
coupled (armchair edge) and symmetric otherwise (zigzag 
edge). The eigenenergies of these potentials have been 
calculated recently. 30 Another recent worlsP3 provided a 
scmiclassical framework to study analytically the edge states 
for both zigzag and armchair edges, but, this approach 
does not furnish the full low energy spectrum. We have 
developped here two semiclassical methods to calculate 



FIG. 21: (a) Red curves: Low energy spectrum for an armchair 
ribbon in a magnetic field in units of t. The width L of the 
ribbon is L — 201ao/2 and the dimensionless magnetic flux / 
is / = 0.00126, so that L/Ib = 9.6. The black curves indicate 
the situations where the classical cyclotron orbits touch an edge, 
(b) Schematic representation of the different situations for the 
position of the cyclotron orbits with respect to the edges of the 
ribbon. 



the energy edge states analytically. A very simple one 
consists in using the Bohr-Sommcrfeld quantization of the 
action related to the effective Schrodinger equations. This 
approximation captures the essential of the physical picture 
except when the cyclotron radius of the effective particule is 
of the order of the distance to the edge. This approach also 
reveals different regions as a function of the energy and the 
distance to the edge, where the skipping orbits are quan- 
tized in different ways. The second method is based on the 
WKB formalism, and accounts for the overlap of the wave 
function when it is close to the edge. Consequently, this 
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more sophisticated analytical approach perfectly describes 
the edge states whatever the distance of the center of the 
cyclotron motion to the edge. In particular, it quantitatively 
describes an interesting region in the armchair case (called 
region T> in the paper) where the energy does not increase 
monotonously. This implies that the drift velocity along 
the edge v2 ri f t = de n /dk y ex de n /dx c may change in sign 
when the position x c varies. These WKB results perfectly 
fit numerical exact calculations. 



Region A : R < \x c \ 



This case, characterized by E < x 2 /2, corresponds to the 



region A in Fig. 22 We write the WKB wave functions in 
both wells and then the matching conditions in x = 0. As 
the potential is symmetric with x, we focus on the left well 
(x < 0). The energy E defines two classical turning points 
at positions x\ — x c — R and xi — x c + R (Fig. 22 1. Within 
the WKB approximation, the connection procedure near X\ 
implies that the wave function in the well reads : 
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Appendix: The WKB method 



(Pl~{x) = J— sm(S(x 1 ,x) + -) 

y/k(X) 4 



(A.3) 



where the left arrow indicates that this wave function 
matches the correct connection procedure near the left turn- 
ing point x\. Ci is a constant and S(xx,x) is the partial 
action in the left well : 



In the whole appendix, we define the classical action as 
half of the one defined in (29). 



In this appendix, we present the detailed calculations on 
the energy levels E n (x c ) within the WKB approximation, for 
the potential V(x) defined as : 



V(x) 



x c f + f 



(x + x c y 



for x < 
for x > 



(A.l) 



The case Vq — has been introduced in section [III B| for the 
study of edge states in zigzag nanoribbons, while the case 
Vo = 1 has been used in section |III C| to describe armchair 
ribbons. 




FIG. 22: The double symmetric well. The two wells are centered 
in ±x c . The two turning points in the left well are x\ and X2- 



S(xi, x) — I dx\l E 



^E-V S ym{x) . (A.4) 



The above expression (A.3) breaks down near the second 
turning point xi — x c + R. Near this point we linearize the 
potential as Vs ym (x) « =^- + R(x — Xi) and the Schrodinger 
equation reads : 



dx 2 



2R(x-xi) tp l (x) = 



(A.5) 



The solution of this equation is given by a combination of 
Airy functions : 

<Pi*(x) = cE;Ai(aa; 4- xq) + J3iQ\{ax + Xo) (A. 6) 
a 3 = 2R (A.7) 

xo=(2R)^ 3 (R-\x c \) (A.8) 

where the right arrow indicates that this wave function must 
obey proper matching conditions near the right turning point 
x\. di and /3i are constants. Inside the well, this wave func- 
tion has the asymptotic expansion : 



<Pi W 



ai . / l , ,3/2 

-TJi Sln I gM 



fn\z\ 
/n\z\ 



i/4 C0S U IZ 



3/2 



with \z\ — —xq — ax = 



k 2 (x) 



(A.9) 



(A.10) 



1. Double symmetric harmonic well 



The argument in the trigonometric functions can be related 
to the partial action between x and xi 



We consider the peculiar case Vo 
symmetric harmonic potential : 

1 



0, that is a double 



v Sym {x) = -(\x\ + x c y 



(A.2) 



For a given energy E. we define the cyclotron radius as 



R = V2E (Fig. 22 1. In addition, in the whole appendix, 



the factor 2 in the definition of the action (29 1 is missing: 
the actions Si and S r are now defined as the half of those 
introduced in Sec. HVl 



\\ Z f 2 = J ' dx \/ E ~ Vsym{x) = Si- S( X1 ,X) (A.ll) 

where Si is the total action in the left well. Therefore ipj*(x) 
can be rewritten in the form : 



<pj*(x) ~ — . cos I S(xt, x) 



V s ' 



+ 7fer i "(*'' :c)+ i~ 5 ') {AA2) 
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with on = otiyf^- (idem for /?j). Then, we impose ip^~(x) 
(pj*{x) inside the well. This implies ai = Ci sin Si and j3i 
Ci cos Si what finally leads to the important relation : 



tanS 1 ; = -g- — 



a i 



, (A.13) 

Pi Pi 

As the two wells are identical, we obtain the same relation 
for the right well : 






(A. 14) 



and, similarly to (A.6), the wave function in the right well 
reads near x — : 

ip^T{x) = a r Ai(— ax + xq) + /3 r Bi(— ax + xo).(A.15) 

The next step is to impose the matching of the two wave 
functions and their derivatives in x = : 

tf*(0) = tf"(0) , d a <f?(0) = d x cptT(0) (A. 16) 

These matching conditions give the two equations : 

diAi(x a ) + /3iBi(x ) = d r Ai(x ) + /3 r Bi(iz; ) 

diAi(x a ) + (3iBi(x ) = -d r Ai'(xo) - /3 r Bi'(x ) 

from which the ratio 2t is extracted as : 

a r _ f (Bi'(so)Ai(xo) + Bi(a;o)Ai'(a:o)) + 2Bi'(x )Bi(x ) 
Tr ~ f2Ai'(x )Ai( a; o)+Ai'(a;o)Bi(.T )+Ai(x )Bi'(xo) 

(A. 17) 
The s ymme try o f the p otential implies Si — S r , so that, 
from (A.13I and iA.UL we have &■ = g- = X. The rela- 
tion (A. 17) becomes a simple polynomial for the unknown 



quantity X : 
Ai'(a;o)Ai(x )A 2 + (Ai'(x )Bi(xo) + Ai(x )Bi'(x )) X 
+ Bi(x )Bi'(x ) = (A.18) 

whose solutions 

Bi'(aro) vA s Bi (^o) 



X c 



Ai'(xo) 



X 



Ai(x ) 



(A.19) 



correspond respectively to the symmetric and antisymmetric 
wave functions. The action Si is quantized by the condition 
tan Si = X S ' AS , so that S t = ir{n + 7) with < 7 < 1, with 
two solutions for 7 : 



l S (E,x c ) 



1 Bi'Qro) 

— arctan -—=- — - 
Tr Ai'(xo) 



-y AS (E,x c ) = 1 -- arctan -4^4 
7T Ai(x ) 



(A.20) 

(A.21) 



In region A, the action is very simpl y related to the energy, 
S = I = nR 2 /2 = ttE, so that from ( |A.2l|A.2l| , we obtain 
the implicit equations : 



E* = n- 



E 



AS 



1 



- — arctan ( — — =- — - 
7T V Ai'(xo) 

1 



arctan 



Bi(x ) \ 
n \ Ai(x )J 

where xq itself depends on the energy : 

1/3 



(A.22) 
(A.23) 



2^2i^ {^2E~-\x c \) 



x ■■ 

These implicit equations are solved numerically and we ob- 
tain the spectrum (region .4) shown in Fig. 23 



b. Region B : \x c \ < R 

This is the region B illustrated in Fig. [22] The expression 
of the WKB wave function matching the correct connection 
procedure near x\ is still given by Eq. (A. 3). The second 



turning point in xi does not exist anymore, and we have to 
know the expression of the wave function near x = 0. It 
reads : 



(pf(x) = a;Ai(y + ax) + /3iBi{y a + ax) (A.24) 

(A.25) 



a 3 = 2 \x r 



Vo 



R 2 



(2M) 



2/3 



(A.26) 



Far from x = 0, the Airy functions can be expanded to 
obtain : 



VT*(?) 






Oil . I I . .3/2 

-^74 sm ( 3 1*1 



|1/4 COS U 10 



3/2 



where we have set : 
M = -2/0 - ax 



k 2 (x) 



OCX 



07 



(2|x c |) 1 / 6 . (A.27) 



The arguments in the trigonometric functions can be related 
to the partial action between x and : 



M 



3/2 _ 



/ dxy/R 2 -x 2 c -2\x c \x + 6 

J x 



Si -S(xi,x) 



where the quantity : 



(R 2 



.2\8/2 



3\Xr 



(A.28) 



(A.29) 



accounts for the step between the potential V(x — 0) and the 
energy E. The wave function <pj*{x) can thus be rewritten 
as : 

<Pl*(x) ~ J— eos(S(x l7 x) + --Si -5) 

+-^=sin(S(xi,x) + j-Si-S). (A.30) 

V k ( x ) 4 

The matching of the two wave functions (pf~(x) — (p~[*(x) 
gives : 

ai 



tan(5; + S) = 



A 



(A.31) 



and a similar expression for the right well. Again we impose 
the continuity of the wave function and its derivative in x = 
and obtain the two equations : 

diAi(y ) + /3;Bi(y ) = d r Ai(y ) + r Bi(y o ) 
diAi'(y ) + (3iBi'(y ) = -d r Ai'(y )-/3 r Bi'(y ) 

from where we extract the ratio a r / ' j3 r . Then, using the 
symmetry 011 — a r and /3; = /3 r , we obtain a polynomial in 
X = ai/ Pi : 

Ai'(y )Ai(y )X 2 + (Ai' (y )Bi(y ) + Ai(y )Bi' (y )) X 

+ Bi(y )Bi'(j/ ) = . (A.32) 
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The only difference with the polynomial (A. 18 1 for the region 
A consists in the substitution xq — > yo 
solutions are : 

Bi'(yo) 



X s = - - 



X AS = - 



Consequently, the 
Bi(ito) 



(A.33) 



Ai'(yo) " Ai(j/ ) 

The action Si is quantized by the condition tanS"; = X S ' AS , 
so that Si = 7r(n + 7) with < 7 < 1, with two solutions for 

7 : 

"1 



7 



S,AS 



(E,x c 



-(arctan X 



S,AS 



S) 



(A.34) 



where a(x) = x— [x], [x] being the next smallest integer. By 
construction < a[x] < 1. 

The last step is to relate the energy to the action. This 
relation is not linear as in region A but it now reads (see Eq. 
32l : 



e e 



1 



sin 20 



(A.35) 



where 9 = arccos x c /R = arccos x c j\JlE has the meaning 
of a angle (Sec. |V| . From this equation, together with the 
quantization condition : 



Si=n{n + j b ' Ab ) (A.36) 

where r y S ' AS (E,x c ) are functions of E and x c through Eqs. 
(A.34 A.26| , we can extract numerically the eigenener gies 

fT 1 '-' onrl Tpy^ OTT/-1 r-\li^f f hnm o c: o fii "ri /"»-H rwi r\+ nr in TTi rv 0*3 



E b and E£ an d plot them as a function of x c in Fig. 



23 



In Fig. ( 24 ) we have also plotted the x c dependence of the 
mismatch index 7 s - (x c ). 
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FIG. 23: Semiclassical spectra E n (x c ) within the semiclassical 
approximation (Dashed lines) 7 = est, (full lines) WKB method. 
The dots are the exact numerical solutions of the Schrodinger 
equation with the potentiel V sy m{x). The two regions are sepa- 
rated by the parabola E = x^/2. For the zigzag problem, we keep 
only antysymmetric states that are the high energy levels (blue 
lines) . 



2. Double harmonic asymmetric well 



In this section we introduce an asymmetry Vq 
potential which becomes : 



1 in the 



rAsym(X) 



2 V "^cj ' 2 

2 \% 1 X c ) 2 



for x < 
for x > 



(A.37) 




FIG. 24: Mismatch index 7 ' (x c ) as a function of x Cl for (blue) 
n = 0, (green) n — 1 and (brown) n = 2. For x c = R, 7 (R) = 
1/3 whereas 7 (R) — 2/3 (horizontal dashed lines). For x c — 0, 
7 S (0) = 1/4 and 7 AS (0) = 3/4. 



and calculate the spectrum E n (x c ) within the WKB approx- 
imation. We now introduce two cyclotron radii Ri and R r . 
The relation between these parameters and the classical en- 
ergy is given by : 



E=^ 



1 



m-i 



(A.38) 



We now have to distinguish three regions delimited by \x c 



Ri and \x c \ = R r and illustrated in Fig. 25 The total action 
in these three regions reads : 



H, r ,ion C ,S = ^Rf , .,.., 



TT 

2 

7T 

2 

R? 



-R„ 



Region D S t = ~Rf + ^ [9 r - ~ sin(20 r )] 



R r 

T 



Rc-ion I) S, = ^i-\ M™l)} + *Y%-\ sin(20 r )] 

(A.39) 



We treat the three regions separately. 




FIG. 25: Double asymmetric well VAsym{x): C, T> and £ refer to 
the three regions discussed in the text. For a given energy E, we 
define two cyclotron radii Ri and R r . 
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3. Region C : R, < R r < \x c 



with 



This region is explicited in Fig. [25] Because of the asym- 
metry of the potential, the action is now different in each 

well : 



Si = -Rf = n(E- 1/2) 



ij r 



-R 2 r = ir{E + 1/2) . 



(A.40) 
(A.41) 



a'i = 2Ri a r = 2R r 

xi = (2i?0 1/3 (i?/-|x c |) 
x r = {2R r ) 1 ' 3 {R r -\x c \) 



(A.44) 
(A.45) 
(A.46) 



Then, by imposing the current conservation and the conti- 
The WKB wave function has the same form as in Eqs. Q), nuit y of these wave functions in x = 0, we obtain the two 
(A. 6 1 and (A. 9 1 written for the symmetric potential so that equations : 



we still have the relations (A. 13 1 and (A. 14 1. Since the po- 



tential is now asymmetric, the expression of the wave func- 
tion near the two inner turning points located in — \x c \ + Ri 
and \x c \ — R r is now : 



a;Ai(x ; ) +/3 ; Bi(x() = a r Ai(x r ) + ,S r Bi(x r ) (A.47) 
diaiAi (xi) + PiaiQi (xi) = —a r a r A\ (x r ) — f3 r a r B\ (x r ) 



(pj*(x) = aiAi(aix + xi) + PiBi(a t x + xi) (A.42) 
ip^T(x)) = a r Ai(—a r x + x r ) + /3;Bi(— a r x + x r )A.43) from where we extract the ratio ~- as : 



,V3 



~ ' Bi'(a; r )Ai(a:,) + Bi(x r )Ai'(x;) (§) J + Bi'(x r )Bi(x;) + Bi(x r )Bi'(x;) (^) 



Br 



21 I Ai 



i'(x r )Ai(x;) + Ai(x r )Ai'(x ; ) (f^) J + Ai'(x r )Bi(x ; ) + Ai(a: r )Bi'(a;j) (f 



1/3 



1/3 



(A. 



Since Si ^ S r , the ratios -= and ~~ are now different. Intro- Then, we insert the relation ( ]A.49 1 into (A. 48) to obtain a 



ducing the total action S t = S r +Si = 2ttE, these coefficients polynomial in X h = ai/fa whose the coefficients are only 
are related as : functions of E and x c : 



a r tan S^ — an 



P r l + ai/fytanSt 



(A.49) 



Hi 
R, 



1/3 



tan. S t ( Bi' (x r )Ai(xi) + [ -^- ) Bi(x r )Ai'(xi) - Ai'(x r )Ai(xi) 

1/3 



R, 



1/3 



Ai'(x;)Ai(x r ) } X? (A.50) 



tanS t Ai (x r )Ai(xi) +Bi' (x r )Bi(xi) + 



Rr 



[Ai' (xi)Ai(x r ) 4- Bi'(x r )Bi(x;)] 



+Bi'(x r )Ai(x,) - Ai' (x r )Bi(xi) + 



Rl 
R, 



1/3 



[Bi(x r )Ai'(xi) - Ai(x r )Bi'{xi)] \ X l 



tan S t Ai'(x r )Bi(x;) + 



1/3 



Ai(x r )Bi'(xO + Bi'(x r )Bi(x,) 



Rl 

R, 



1/3 



Bi(x r )Bi'(x/) = 



with St — 2nE. Because of the asymmetry of the potential, 
the coefficients of the polynomial are much more complicated 
than for the symmetric case (Eq. A. 18). We check that for 
Rr = Ri and xi = x r = xo we recover the symmetric case 
(Eq. A. 18 1. The polynomial (A.50) has still two solutions 



X, and X^ s , from which we obtain the action Sr. It has 



v i > 



the from Si — ir(n + 7) with : 



1 

7T 



"ff(E,x c ) — — axctanXf 



Ji(E,x c ) = - arctanX 

7T 



AS 



(A.51) 
(A.52) 
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Since the energy E is simply related to the action E = 1/2 - 
Si/ iv, we finally obtain the two implicit equations : 



E* = n 

E AS 



n 



i + iarctanJff (A.53) 

§ - ^arctanA/ 45 (A.54) 



from which we obtain the energy levels E n (x c ) in the region 
Ri < R r < \x c \. 



The next step consists in the matching of the Airy functions 
with the WKB approximation valuable inside a well. The 
treatment for the left well has been performed in the para- 
graph [Ta] of this appendix and gives: 



, c Oil (Ml 

tan bi = — = -=- 

Pi A 



(A.61) 



4. Region V : R t < \x c \ < R r 

The region T> represented in Fig. [25] does not exist for 
the symmetric potential Vs ym (x). We have now to linearize 
the potential around — \x c \ + Ri for the left well and around 
x = for the right well. The solution inside the left well 
takes the familiar form : 

¥(x) = diA\(aix + xi) + B~iQ\(aix + xi) (A.55) 



a? = 2Ri 

xi = (2Ri) 1 / 3 (Ri-\x c \) 

whereas the one inside the right well is 



(A.56) 

(A.57) 



ty r (x) = a r Ai(y r — a r x) + j3 r Bi(y r — a r x) (A. 58) 
a P = \2x c \ (A.59) 



Mr 



\2x, 



1 2/3 



(A.60) 



whereas the calculation for the right well has been made in 
the paragraph 1 1 b| and gives : 



(A.62) 
(A.63) 



tan(5V + § r ) 
where S r 



a r 


a r 

J r 


(R 2 r - 


2^ 3 / 2 



3 b, 



The current conservation and the continuity of the wave 
function at x = gives the two equations : 



aiM(xi) + ABi(arj) = a r Ai(y r ) + (3 r Bi(y r ) 
diaiA\(xi) + (3iaiBi'(xi) = ~d r a r Ai' \y r ) - {3 r a r Bi'(y r ) 



from where we extract the ratio %=•: 



- r f \Bi'(y r )Ai(xi) + (g|) Bi( Vr )Ai'(a:i)J + Bi' (y r )Bi( Xl ) + (^) ' Bi(y r )Bi> (xi) 



& f (Ai'(y r )Ai(^) ' I ' 



1/3 



1/3 



.i^-.v ,- v — j , Ai(y r )Ai / (a:i)j+Ai / (v P )Bi(aJi) + ( 1 ^|-) Ai(y r )Bi'(^) 

Introducing the total action St = Si + S r and using ( JA.61 1 and (A.62 1 we obtain the relation : 

a r _ tan(St +Sr)-% _ v - Xi 
Tr ~ l + tan(St + 6 r )% ~ 1 + vXC 



(A.64) 



(A.65) 



where we have introduced v = tan (S t + S r ). We inject this last relation into (A.64 1 to obtain again a polynomial in Xi : 



v Bi'fo r )Ai(a:j) 



Ri 



1/3 



Bi(y r )Ai'(xi) ) - Ai'(p r )Ai(x,) - ( ^ ) Ai' (xi)Ai(y r ) } Xf 



Ri 



1/3 



(A.66) 



v [ Ai'(y r )Ai( Xl ) +Bi'(y r )Bi(x l ) 



+Bi'{y r )Ai(xi) - Ai'(y r )Bi(xi) + 



1/3 



Ri_ 

\x c \ 



1/3 



[Ai'{xi)Ai{y r )+Bi'{xi)Bi{y r )] 



Ri 



1/3 



[Bi(y P )Ai'(a:,) - Ai(i/ P )Bi'(xi)] \ *i 



1/3 



+v I Ai'(y r )Bi(xi) + [ A ) Ai(y r )Bi'(xi) ) + Bi' (y r )Bi( Xl ) + [A ) Bi(y P )Bi'(^) = 0. 



Note that this polynomial satisfies the continuity of the spec- 
trum on the parabola for which x c — R r . Indeed, in this case 
we have 8 r = and x r = y r — so that the two polynomials 



(A. 50 1 in region C and (A.66 1 in region V coincide and give 



the same solutions. From here, the work is exactly the same 
than for the lower region. Therefore, we have just to take 
the two solutions Xf and Xf- S of this polynomial and insert 
them into (A. 51 1 to obtain the left MI in the intermediate 
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region. Then we insert them into (A. 53) and (A. 54 1 to ob- 



tain the implicit equations in E n and x c from where we can 
extract the spectrum E n (x c ) in the region T>, characterized 
by the relation R[ < \x c \ < R r . 



5. Region £ : \x c \ < Ri < R r 

The last region £ satisfies \x c \ < Ri < R r and is illustrated 
in Fig. |25j As for the high energy region of the double sym- 
metric well, we decompose the well into two parts from x = 0, 
and keep working with Si that is the action in the well where 
x < 0. The strategy is still the same, and the calculations 
are very similar to those in section 1 1 b| We express the ra- 
tio a//3 from two ways. First, we obtain a relation of the 
type (A. 31 ) thanks to the matching of the WKB wave func- 



tion with the solutions of the linearised Schrodinger equation 
around x — that we express as : 

¥(x) = diAi(aix + yi) + PiB\(aix + y t ) for x < 
^ r (x) = d r A\(y r — a r x) + f3 r Bi(y r — a r x) for x > 



Ui 



•*-> f. -*W *** r> ■*«»■ 



1 2s, 



,2/3 



Vr 



\2x, 



,2/3 



a? = al = \2x r 



The difference here comes from the asymmetry Vo = 1 so 
that the relation (A. 31 1 becomes : 



tim(Si + 8i)=X l 

tan (S r + 5 r ) = X 

s r = (^-^ 

3 \xc. 



3/2 



a r 

~~Jr 



(R? 



,3/2 



3 \Xr. 



(A.67) 
(A.68) 

(A.69) 



In order to obtain a polynomial in terms of Xi in the high 
energy region, we use the lattest relations to express X r as 
a function of X\ : 



J\ r — 



J\ r — 



tan (St + S r ) — tan Si 
1 + tan Si tan (St + 5 r ) 

(v tan 5i — 1) %- + v + tan Si 
(v + tan Si) ^ - (v tan Si - 1) 

u Xl+T 
T Xi-O 



(A.70) 



where we have introduced a = v tan Si — 1 and r = v-\-tanSi- 
Then, we use the continuity conditions at x = to obtain : 



yi. ,. 



o, = Xi [Bi'(y r )Ai(yi)+Bi(y r )Ai'( yi )] + Bi^y^Bi^) + Bi(y r )Bi , (y t ) 
p r = Xi [Ai'(y r )Ai( yi )+Ai(y r )A\!( yi )] + Ai' (y r )Bi( yi ) + Ai(y r )Bi' ( Vl ) 

>From Eqs. (A.70) and ( |A.71[ ), we obtain a polynomial for Xi : 

{r (Bi'(y r )Ai(yi) + B% r )Ai'(yO) + a (Ai'(y r )Ai(yO + Ai'( W )Ai(|fc.)) } X? 
+ {r (Ai'(y r )Ai( W ) + Bi'(y r )Bi(y ; ) + Ai'(yOAi(y r ) + Bi'( W )Bi(ifr.)) 
-a (Bi 1 (y r )Ai(yi) - Ai'(y r )Bi(j/,) + Bi(y r )Ai'(y ; ) - Ai(j/ r )Bi'(j/ z ))} X t 
+t (Ai'(y r )Bi(yi) + Ai(y P )Bi'(y,)) - o (Bi'(y r )Bi( W ) + Bi(y r )Bi'(yO) = 



(A.71) 



(A.72) 



v. 



Again, it is easy to check that for |x c | = Ri, we have r 

a = — 1 and Xi = yi — so that the polynomials ( 1 b I 



in region V and (A.72 1 in region £ coincide, what assures 



the continuity of the spectrum between these two regions. 
To obtain the implicit equations and extract the spectrum 
E n (xr) in this region, we inject the solutions X S I AS into 
(A.67 1 to obtain an expression of the mismatch indexs j s / AS 
like in region B (see Eq. (A.34) with 5 —> Si) and write 



explicitly the action Si as a function of the energy and x c as 
in the expression ( 32 1 with R—}Ri. 
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